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Abstract 

Energy distributions of high frequency linear wave fields are often modelled in terms of flow or transport 
equations with ray dynamics given by a Hamiltonian vector field in phase space. Applications arise in under- 
water and room acoustics, vibro-acoustics, seismology, electromagnetics, and quantum mechanics. Related 
flow problems based on general conservation laws are used, for example, in weather forecasting or molecular 
dynamics simulations. Solutions to these flow equations are often large scale, complex and high-dimensional, 
leading to formidable challenges for numerical approximation methods. This paper presents an efficient and 
widely applicable method, called discrete flow mapping, for solving such problems on triangulated surfaces. 
An application in structural dynamics - determining the vibro-acoustic response of a cast aluminium car body 
component - is presented. 



1 Introduction 



Wave energy transport in the high frequency limit can be described in terms of the underlying ray dynamics, 
neglecting interference and other wave effects. The problem is thus reduced to tracking ray densities in 
phase space and becomes part of a wider class of mass, particle or energy transport problems driven by an 



underlying deterministic velocity field. Applications are found, for example, in fluid dynamics ( Celani et al 



( 2004 1), weather forecasting (Sommer & Reich ( 2010 1) , linear wave dynamics or in general in describing the 



evolution of phase space densities of a dynamical system. The governing equations for the particle or ray 
densities can often be written in terms of conservation laws (LeVeque (19921), an example is the Liouville 
equation describing the evolution of ray densities in phase space for Hamiltonian flows. This class of problems 
is particularly interesting in the context of approximating linear wave equations in terms of their underlying 
ray dynamics in the short wavelength limit 

Numerical approaches for solving transport problems of this kind are typically based on the method of 
characteristics, that is, the solutions are found along trajectories or rays determined by the underlying vector 
field. The flow equations can be formulated in terms of a linear propagator, the so-called Frobenius-Perron 



(FP) operator (see, for example, Cvitanovic et al. ( 2012 1) , which describes the evolution of phase space 
densities in time. Numerically efficient methods for solving flow problems in more than one dimension for a 
wide class of physically relevant systems are still non-existent. A variety of techniques have been developed 
based on a FP-operator approach, however, all with a fairly limited range of applicability. Difficulties arise 
due to the high-dimensionality of the phase space and the singular nature of the operator describing the 
underlying deterministic dynamics. One approach for dealing with such problems is Ulam's method (see 
e.g. Ding & Zhou (1996)), which is based on subdividing the phase space into distinct cells and considering 



transition rates between these phase space regions. Other methods include wavelet and spectral methods for 



the infinitesimal FP-operator (Junge & Koltai (20091; Froyland et al. (20131) and periodic orbit expansion 



techniques ( Lippolis & Cvitanovic ( 2010 1; Cvitanovic et al. ( 2012 \). The modelling of many-particle dynamics 



such as protein folding, has been approached using short trajectories of the full, high-dimensional molecular 



dynamics simulation to construct reduced Markov models ( Noe et al. ( 2009 1). For a discussion of convergence 
properties of the Ulam method in one and several dimensions, see Bose & Murray (20011 and |Blank et al] 
( j2002[ ), respectively. 

More direct methods are based on tr acking swarms of trajectories in phase space often referred to as ray 
tracing, see for example Cerveny (20011. Methods related to ray tracing but tracking the time-dynamics of 



interfaces in phase space, such as moment methods and level set methods, have been developed by |Osher fc| 



Figure 1: (a) Discrete flow mapping on a mesh of a thin aluminium shell (shock tower) from a Range Rover; (b) 
the boundary flow map (pij between a pair of adjacent triangles. (Online version in colour.) 



Fedkiw (20011, Engquist & Runborg (20031, Ying & Candes (20061 and Boon et al. (20071 amongst others. 



They find applications in acoustics, seismology and computer imaging, albeit restricted to problems with 
few reflections; for an excellent overview, see |Runborg| ([2007 ). In the following we focus on ray-tracing 
approximations of linear wave problems, although the methodology developed here can be used in a more 
general context. 

Ray tracing and tracking methods often become inefficient when considering stationary (wave) problems in 
bounded domains, or in general for ray tracing problems including multiple scattering trajectories and chaotic 
dynamics. An example is the wave field in a finite cavity driven by a continuous monochromatic excitation. 
Here, multiple reflections of the rays and complicated folding patterns of the associated level-surfaces often 
lead to an exponential increase in the number of branches that need to be considered. It is thus necessary 
to employ approximation methods for the associated ray-tracing solutions, often based on ergodicity and 
mixing assumptions of the underlying ray dynamics. A popular tool amongst the mechanical engineering 
community in the context of vibro-acoustic modelling is statistical energy analysis (SEA) (see for example 



Lyon (19691, Keane (19921 and Lyon & DeJong ( 1995[ )). A related method for electromagnetic fields is the 



random coupling model, which makes use of random field assumptions (see Hemmady et al. (2012 i). In SEA 
the structure is subdivided into a set of subsystems and ergodicity of the underlying ray dynamics as well 
as quasi-equilibrium conditions are postulated. The result is that the density in each subsystem is taken to 
be approximately constant leading to greatly simplified equations based only on coupling constants between 
subsystems. The disadvantage of these methods is that the underlying assumptions are often hard to verify 
a priori or are only justified when an additional averaging over 'equivalent' subsystems is considered. The 
shortcomings of SEA have been addressed by Langley ( 1992 I, Langley & Bercin ( 1994 1 and more recently in 



a series of papers by Le Bot ( 1998 2002 20061 



A computational method called dynamical energy analysis (DEA), which is based on an FP-operator 
approach, has been introduced by Tanner (20091 and further developed in Chappell et al. (2011 1. DEA is an 



Ulam-type method subdividing configuration space into smaller subsystems, but using a basis approximation 
(spectral method) to obtain an improved resolution of the phase space density. By increasing the resolution in 
both the position and momentum variable, one systematically interpolates between SEA and full ray tracing 
thus relaxing the underlying ergodicity and quasi-cquilibrium assumptions in SEA. A more computationally 
efficient approach using a boundary element method for the spatial approximation has been applied to both 
two and three dimensional problems in Chappell et al. (20121 and Chappell & Tanner (20131. A major 



advantage of DEA is that by removing the SEA requirements of diffusive wave fields (equivalent to the 
ergodicity assumption) and quasi-equilibrium conditions, the choice of subsystem division is no longer critical. 
The resulting increase in flexibility in this choice leads to a much more widely applicable method. 

In this work, we exploit the freedom in choosing the subsystems by extending DEA towards solving 
the Liouville equation (or other hyperbolic equations) on meshes. Meshed surfaces are replete in numerical 
simulation problems thanks largely to the huge popularity of finite element methods. Considering the elements 
of a mesh as our basic domains therefore automatically renders our techniques applicable to a wide class of 
problems including the ability to handle complex geometries; an example of which is given in Fig. [TJl A 
highly efficient solution procedure can be developed for triangulated surfaces since the problem is simplified 
to considering the local ray dynamics in flat planar regions. In addition, solving phase space flow equations 
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on meshes makes it possible to consider geodesic dynamics on curved surfaces or ray dynamics in non- 
homogeneous media. This opens up the range of possible applications enormously, for example, in modelling 
high frequency vibrations of thin elastic shells (see Norris (19951 and Yang et al. (1996 1 ) or in underwater 
acoustics (Jensen et al. ( 1993 1) . 

The main outcome of this paper will be the introduction and development of the Discrete Flow Mapping 
(DFM) technique, a new and efficient method for solving stationary phase space flow equations in complex 
domains. DFM has similarities to both finite volume methods and boundary integral methods, combining 
the advantages of each. One achieves the reduction in dimensionality to the boundary of each sub-domain 
characteristic of boundary integral methods, which in phase space is a reduction by two. In addition one can 
treat non-homogeneous domains as in finite element and finite volume methods by applying an approximation 
of local homogeneity and refining the mesh to improve this approximation. A great advantage of DFM is 
that it can be applied directly to existing finite element models with relatively coarse meshes, and provides 
a solution for the high frequency case which automatically contains the geometric details absent from SEA- 
type methods. In addition, the flow directivity can be resolved, unlike in an SEA model, where ergodicity is 
assumed. 

The paper is structured as follows: we will give an integral equation formulation for the stationary Liouvillc 
equation on triangulated surfaces in Sec. [2] We then detail the implementation of DFM in Sec. [3] using the 
linearity of the integral operator, we approximate the solution using a mixture of boundary element and 
spectral methods. We also discuss the treatment of reflection and transmission on interfaces, such as at edges 
or due to abrupt changes in material parameters. Finally, in Sec.[4]DFM is applied to model the phase space 
flow on a sphere and the geodesic flow on an irregularly shaped car body part, demonstrating its power and 
efficiency in practice. 



2 Integral Equation Formulation 



We will focus on triangulated surfaces Q = U -=i ^ consisting of Nq triangles Qj, j — 1, ...,Nn such 

as depicted in Fig. [Ik. We consider piecewise constant Hamiltonians of the form Hj(r,p) — Cj\p\ = 1 in Qj 
describing the energy of a flow, where c, is the flow velocity for r £ Qj and the momentum coordinate p lies 
on a circle of radius c~ l . This Hamiltonian is associated to the Helmholtz equation with inhomogeneous wave 



velocity c(r) (see Runborg (2007)). In Sec. 3) |3.2 l we also consider vectorial wave equations, such as plate 
equations, which include different wave modes. In this case the wave propagation needs to be characterised 
by more than one Hamiltonian per triangle j, that is, we need to consider Hamiltonians H l j — c l j\p\, where 
I refers to the mode type. We restrict our discussion to the scalar case in the following for simplicity of 
notation. 

Let us denote the phase space on the boundary of the triangle Qj as Qj = dQj x (—cj 1 , cj 1 ). The 
associated coordinates are Xj = (sj,pj) £ Qj with Sj parameterising dQj, the boundary of the jih triangle, 



and pj e (— c ■ 



parameterising the component of the inward momentum (or slowness) vector tangential 



to dQj. We denote the boundary flow map as ipij : Qj — > Qi which takes a vector in Qj and maps it under 
the flow given through Hj to a vector in Qi, see Fig.[lj3. Note that tfiij is generally only defined on a subset of 
Qj, namely the preimage of Qi, that is tpT^ (Qi) C Qj. This preimage is empty if Qj and Qi are not adjacent 
(see Fig. [l}. 

The stationary density p(Xi) on Qi, i = 1, Nq, due to an initial boundary distribution p 1 - ' on Qj, 
j = 1, Nn, may be computed using the following boundary integral equation (see Tanner (2009), Chappell 



et al. (2012) and Chappell fc Tanner ( 2013 1), 



(I-B)p(X i ) = pV\X i ) 



where 



Bp(Xi) := Y, 



Kr(Xi,Xj)p(Xj)dsjdpj 



(2) 



and Kr describes the propagation of the flow. Here we consider purely deterministic flows with Kr(Xi, Xj) = 
w(Xi) S(Xi — ipij(Xj)). A diffusion component may be added by replacing the ^-distribution with a finite- 
width kernel. In general, the weight function w contains reflection/transmission probabilities at boundaries 
and a dissipative term of the form exp(— pL), where L is the length of the flow trajectory and p is a damping 
coefficient. For the case w = 1, the transfer operator B is of Frobenius- Perron type with a maximum 
eigenvalue one. In this case the problem, Eq. |l]), becomes ill-posed and thus we consider the case p > 
henceforth. One can also obtain a well-posed problem using other forms of dissipation, such as an absorbing 

( 2012[ )). Note that since the flow map only maps to neighbouring 



or open boundary region (Chappell et al 



triangles, a matrix representation of B over the whole of UiJl Qi ^ s m general sparse. Equation |l| is a 
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Figure 2: Admissible ranges (s TO ,- n , s max ) for Sj when a ray travels from edge I of triangle j with fixed direction 
coordinate pj = sm(9)/cj to an opposite edge a: (a) Sj can take all values in (0, L) on the edge I (i.e. s m i n = 0, 
Smra = L)\ (b) e (0, s max ); (c) e (s min .L). (Online version in colour.) 



boundary integral formulation for the stationary Liouville equation as d etailed in Chappcll & Tanner ( 2013 1 
The derivation of p' ' for a high frequency point source is also given in Chappell & Tanner (20131. 



3 Implementation of discrete flow mapping 



3.1 Discretisation 

Typically it is assumed that f2 consists of a large number of triangles iVn describing a geometrically complex 
domain. A discrete approximation of p is sought in phase space coordinates on the triangle boundaries. The 
spatial approximation is given by taking polynomial basis functions on each triangle edge a with a = 1, 2, 3. 



That is, we split the approximation at corners, see Chappell et al. (2011 1 for further discussion on the benefits 
of this choice. In the following, we use piecewise constant functions on the triangle edges for simplicity. In 

-1 ) only. It 



contrast to the position coordinate, the momentum coordinate has support on the interval (— c~ , 
is therefore proposed to employ a Legendre polynomial basis approximation in this coordinate (see Chappell 
et al. ( 2012 1) . A key advantage of these choices is that the integrand in the operator |2j remains a very 
simple function of the position argument and the corresponding integral can be performed analytically. This 
dramatically reduces the costs of evaluating |2j| compared to the implementation used in |Chappell et al. 
( 20121 ). 

The overall approximation on Qi for i = 1, Nn is then of the form 



3 X p 

p T {Xi) w ^2^2p(i, a ,i3)b a (si)P[i(pi) 



(3) 



where N p is the order of the momentum basis expansion. The momentum basis functions are given by 

Pfs(j?i) =V^Pp(dPi), (4) 

where Pp is the Legendre polynomial of order /3. The piecewise constant spatial basis functions are given by 
b a {si) = 2~ 1//2 /V A a for Si on the edge a, and zero elsewhere. Here A a is the length of the edge a £ {1,2,3} 
of Imposing a weak form of the integral equation |l]) using the orthonormal inner product for Legendre 
polynomials < • , • > yields 



where / is the identity matrix, 



P(i,a,P) 



(/-B)p = p u 



p {0) {s l ,p i ),b a (s i )Pp(p i )) , 



(5) 



(6) 



-B(i, Q ,/3),(j,i,m) = (B (bi(si)P m {pi)^ ,b a (si)Pp(pi)J 
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Figure 3: T-joint as an example of a junction between 3 plates. The dashed line indicates an incoming wave, 
the solid lines represent outgoing bending, shear and pressure waves. (Online version in colour.) 



and the vectors p, p° have entries given by P(i :a ,0)> P% a $)■> respectively. Expanding the inner product in (pi 
and using the definition of the transfer operator B pi, the discretised transfer operator B acting on Qi for 
alii = 1, ...,Nn may be expressed as 

i,a,{3) ,(j,l,m) (7) 

_ 2m + 1 / / p fj (p.}i )a ( s ^K T (x t ,Xj)Prn(pj)bi(sj)dsjdpjdsidpi 
4 JQi JQj 

= 2m 4 + 1 / W^jj^JJP^^fX^ila^^JJPmtpjJ^SjJdSjdpj. 

Here we write ifiij = ((pfj, vf, ) to denote the splitting of the position and momentum parts of the boundary 
map. Using the properties of the weight function and the spatial basis functions, Eq. Q simplifies to 

B(i :a ,0),(i,l : m) = (8) 

""' ~ 1 r m "" &, £Ww))e-^ iC ' ,,v « ( ^ )) ^(^(^i))d»jdpj- 



l (p j ,a,l) 



Here, w rt is equal to the transmission probability w t when i ^ j, and is equal to the reflection probability 
w r = 1 — Wt otherwise. Higher dimensional transmission/reflection matrices arise in the case of multi-mode 
dynamics as discussed in the next section. Also Ly is the length of the trajectory from the point represented 
by Sj to (plJXj), and (s m i n , s max ) is the admissible range of values for Sj. Restriction to this range is 
necessary to ensure that a ray starting on edge I of triangle j with direction coordinate Pj will be transferred 
to a particular adjacent edge a of triangle i as shown in Fig. [2] Note that for flat polygonal elements such 
as triangles, <Pij(Xj) in fact only depends on pj, and hence only the damping term in equation (|8]) retains 
dependence on Sj . The inner integral thus has a simple form and can be computed analytically, leaving only a 
single integral to be evaluated numerically. It is this step which leads to vast improvements in the computing 
times, and enables us to consider many thousands of elements with high order approximations in direction 
space. Two key points are: 

• Triangulation is used to ensure that the elements of a piecewise linear mesh are flat. This enables us to 
deal with complex enclosures/surfaces whilst keeping the system locally simple. 

• Due to the semi-analytic phase space integration, working with many thousands of subsystems and high 
order direction space approximations becomes tractable. 

It is worth pointing out that the analytic integration method is not restricted to triangular elements, but 
works for any flat polygons. The general idea of the DFM is related to the work in |Chappell et aE| ( |2011[ ) and 
Chappell fc Tanner] |2013 ) . However, in these papers, the DEA method was developed for general subsystems 



of arbitrary shape and the double integrals in the equations equivalent to Eq. (|8| needed to be computed 
entirely numerically. Computing the matrix B is then time consuming and the computations were restricted 
to a small number of subsystems (typically 10 or less). 

3.2 Transmission at interfaces and ray tracing on curved surfaces 

For modelling high frequency wave propagation it will be necessary to take into account transmission (and 
consequently reflection) at interfaces with abruptly changing material parameters or due to edges and branch 
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lines. The latter will arise, for example, when plates are connected along a common junction such as depicted 
in Fig. [3] These effects are included in our approach through transmission/reflection probabilities at triangle 
boundaries coinciding with the interface. The reflection/transmission coefficients are obtained from local 
wave solutions at the interface, incorporating the dependence on the direction of the incoming and scattered 
waves. A number of scenarios important in the context of both wave propagation and ray tracing are outlined 
below. 

Scalar wave equations: For scalar wave equations and simple planar interfaces, the transmission coef- 
ficients are derived from the wave equation using continuity of the wave function and its normal derivative 



for an incident plane wave (see DeSanto (1992)). One thus obtains 



w t (ki, kj 



4(ki/kj) cos(#j) cos(#i) 
((ki/kj) cos^-) +cos(6»0) 2 



(9) 



Here, kj and ki are the wave numbers in the elements containing the incoming and outgoing rays, respectively, 
with ki = ui/ci and uj is the angular frequency. A change of wave vector at the interface may be due to a 
change in the material parameters for example. The angle of incidence 9j of the incoming ray with respect to 
the normal at the element boundary is simply arcsin(cjp_, ). The direction of the outgoing ray 8i is determined 
using SnelFs law. 



Elastic waves in coupled plates: A more general case is represented by interfaces forming junctions 
between plates of varying thickness and meeting at arbitrary angles. The transmission/reflection coefficients 
for a set of plates coupled at a common interface (such as depicted in Fig. [3| may be computed using 



the methods presented in Langley & Heron (19901 and Craik et al. (20041. In particular, we consider the 



connection between plates as line junctions, that is, the interior properties of the junction are not modelled 
and the mass and moment of inertia are neglected. Let us consider a line junction which couples n different 
(for simplicity) semi-infinite plates. The boundary conditions at the line junction correspond to dynamic 
conditions involving stresses, moments and kinematic conditions for the displacement and rotation of the n 
plates. To construct the transmission coefficients, we calculate the response of the system with respect to 
excitation by an incoming plane wave. The incoming wave has a fixed wavenumber and a characteristic mode, 
that is, it is of bending (b), pressure (p) or shear (s) type. The outgoing waves typically have components in all 
n plates and are a mixture of all mode types. Evanescent modes may be included to complete the description. 
Possible material differences between the plates can lead to different wavenumbers in different plates. For 
a given forcing with a particular incoming mode in a particular plate, we can solve for the unknown modal 
coefficients in all plates. In practice, we find the transmission probabilities directly by calculating the ratio 



of outgoing to incoming normal power fluxes. A detailed description can be found in Chappell et al. (20131 



Curved surfaces: In the case of a geodesic flow on a smooth curved surface it is necessary to mimic 
the geodesic paths on the corresponding triangulated surface. We employ the theory outlined in |Kimmel fc| 
Sethian ( 1998 1 and Martinez et q/.| ( |2005[ ), making use of the fact that for ray paths not intersecting vertices on 



triangulated surfaces, the notions of shortest and straightest (discrete) geodesic arc equivalent (see Martinez 
\et fflZ.|p005[ )). Hence the straightest geodesic choice of 8j — 8i will approximate the direction of the geodesic 
flow on homogeneous regions of the triangulated surface, and abrupt changes of surface thickness or material 



will result in a Snell Law effect (see Kimmel & Sethian ( 1998 1). A suitable choice of quadrature method here 
ensures that the rays never pass through vertices of the triangulation (although a point source may lie on a 
vertex), for example Gaussian quadrature rules where endpoints are never used as abscissae. 



Curved shells - curvature corrections: For modelling the vibration of thin shells we need to consider 
ray tracing on curved surfaces, where the dynamics are derived from thin shell theory (Norris (1995)). In a 
lowest order approximation, curved rays again follow the geodesies of the surface ( Yang et aZ.| ( |1996| ); |Tanner 
& S0ndergaard (20071). This approximation derives from applying thin shell theory in the high frequency 



limit as in Norris & Rebinsky ( 1994 \ and Norris ( 1995 1, which is valid for wavelengths shorter than the radii 



of curvature of the shell, but larger than the thickness. Corrections to the geodesic ray approximations need 
to be considered if the local radii of curvature are of the same order as the wavelength. It is possible to 
construct modified ray paths from the dispersion curves given by thin shell theory (see |Norris fc Rebinsky 



( 1994 1), but this requires a detailed knowledge of the local curvature. In the interests of keeping the model as 
simple as possible we will follow a different approach here. We treat the meshed structure as a set of plate-like 
elements and estimate reflection/transmission properties due to the finite angle between mesh elements using 
the plate-junction theory sketched above and detailed in |Langley fc Heron (19901 and Craik et al. (20041. 



This enables us to mimic wave barriers due to regions of high curvature as demonstrated in the next section. 
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Figure 4: Ray density (logarithmic scale) on a sphere computed using discrete flow mapping showing the source 
point (upper) and antipode (lower) for different orders of direction space basis (N p ) on a triangulated sphere with 
5120 triangles. The plots on the right show a representation of the exact ray density over the same triangulation. 
(Online version in colour.) 



In each of the cases considered above, reflection/transmission coefficients are incorporated as part of the 
weight function w in the boundary integral kernel Kr in Eq. Q, or its finite dimensional approximation. 
We assume that the transmission coefficients depend only on the incoming momentum p s via the angle of 
incidence; the outgoing angle is given by Snell's law taking into account refraction due to differences in the 
wave speed across the interface. 



4 Numerical examples 

We will demonstrate the efficiency and flexibility of DFM with the help of two examples. Firstly, we determine 
the ray density produced by a point source on a sphere, where an analytic solution is available for verification. 
Secondly, we compute the energy density distribution on the curved surface of a cast aluminium Range Rover 
body part and compare with solutions of the corresponding wave equation obtained using the finite element 
method (FEM). 



4.1 Flow on a sphere 

The ray density p generated by a geodesic flow emanating continuously from a point source on a sphere can 
be determined analytically. As a function of the polar angle (j>, with source point = 0, the ray density is 
given by: 

pW= (l- e -^)sinW (10) 
Here, \jl is the damping coefficient introduced in Sec. [2] and C is a constant depending of the strength 



of the source. The derivation of equation ( 10 1 follows simply from the fact that the ray density on the 
sphere is inversely proportional to the element of surface area. The exponential terms result from damping 
contributions of the form exp(— /j,L); summing over the trajectory lengths L at any point on the sphere results 
in a geometric series with a contribution each time the ray orbits the sphere. 

The sphere is thus a good candidate for verifying our approach on an approximately spherical triangulated 
surface. The dynamics on the sphere are integrable and the exact solution for the density contains a singularity 
at both poles (<j> = and ir) as the rays are unidirectional along great circles passing through the source 
point. This example is therefore slightly atypical since ergodic or mixing ray dynamics in complex geometries 
generally lead to more smoothly distributed ray densities. The sphere is thus a true challenge for DFM, which 
due to the finite order basis approximation will always incorporate diffusive behaviour and consequently a 
smoothing of singularities. 

Fig. [4] demonstrates that DFM can deal even with such a singular case. Higher order implementations of 
the basis approximation in direction space lead to an improved resolution of the refocussed singularity. Table 
[I] gives the mean relative error averaged over the upper hemisphere containing the source point shown in the 
upper row of Fig. [4] The results are given for Delaunay triangulations with differing numbers of elements 
and for different orders of direction space basis. The results are computed at the centroid of each triangle 
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Table 1: Mean relative errors computed at triangle centroids and averaged over the upper hemisphere of a 
Delaunay triangulated sphere with Nq triangles and a direction space Legendre polynomial basis approximation 
of order N v . 



No 


N p 


Mean Relative Error 


320 


4 


0.1606 


320 


6 


0.1129 


320 


8 


0.1142 


1280 


8 


0.08704 


1280 


10 


0.08884 


5120 


10 


0.06648 


5120 


12 


0.06212 


5120 


14 


0.06100 


20480 


14 


0.05116 


20480 


16 


0.05012 



1.95 

1.9 
1.85 

1.8 
1.75 

1.7 
1.65 

1.6 
1.55- 

1.5 
1.45 




Figure 5: Kinetic energy density on a thin aluminium shell estimated using an averaged full wave finite element 
model (left) and discrete flow mapping (right) for 3% hysteretic damping. (Online version in colour.) 



and compared against the exact solution for the same value of <j>- For the results presented here we have 
taken // = 1 and C = (800tt 2 ) - , which is the scaling for a unit excitation of the Helmholtz equation with 
k — 1007T. The factor is derived from the high frequency asymptotics of the 2D fundamental solution for the 



Helmholtz equation (see Chappell & Tanner ( 2013 1) and matching the asymptotics with equation ( 10 1 as the 
distance from the source becomes small. 



4.2 An Application in Vibro-acoustics 

In this section we consider the transport of high frequency flexural or bending wave energy through curved 
structures using thin shell and high frequency ray models. At present there is wide interest in modelling 
aluminium shells in the automotive industry due to a desire for lighter, and hence lower emission vehicles. In 
particular, large molded aluminium components are replacing more traditional multi-component beam-plate 
constructions, which has the additional advantage of eliminating problems due to fatigue at joints. High 
frequency vibro-acoustic models based on an SEA treatment will be unsuitable in these circumstances, since 
complex geometrical features are not included in SEA. In addition, a subdivision of the model into subsystems 
is not clear cut for such castings as all components are well connected, see for example the structure in Fig. 
[T^i representing an aluminium shock tower from a Range Rover. 

Discrete flow mapping can overcome these problems, since it can be easily applied in the framework of 
existing grids for finite element models, requires no choice of subsystem division and incorporates the full 
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Figure 6: Kinetic energy (KE) response against hysteretic damping level. The KE is averaged over the blue 
circular receiver regions. DFM: solid blue line with 'o' markers; FEM: dotted red lines with '•' markers; mean 
FEM: solid green line with 'x' markers. (Online version in colour.) 



geometry and directionality of the energy flow. The response of a thin molded aluminium car component 
(Range Rover shock tower) to a point force applied normally to the surface is modelled and compared with 
a FEM approximation for the full wave model performed using Nastran. Figure [l^, shows the problem setup 
including the mesh and indicating the local DFM map <pij . In order to maintain a tractable model size for the 
finite element calculation and to study frequency ranges of industrial interest, the computation is performed 
at frequencies between 8 kHz and 10 kHz. This approximately corresponds to a third of an octave band 
centred at 9 kHz. A key assumption of the thin shell theory is that the radius of curvature is large compared 
to the wavelength (see |Norris fc Rebinsky| ( |1994[ )), which is only partly satisfied for this structure in the 
frequency band considered. The wavelength is typically around 3 to 6 cm depending on the shell thickness. 
We therefore employ a modified geodesic ray tracing technique by incorporating reflection/transmission due 
to finite angles between mesh elements as described in Sec. |3 pT2| ). This leads to geodesic trajectories without 
reflections in relatively flat areas, but incorporates wave barrier behaviour in regions of high curvature. 

Fig. [5] shows the DFM result (bending mode only) compared with the Nastran solution for the kinetic 
energy averaged over 41 evenly spaced frequencies spanning the prescribed range and with a hysteretic 
damping level of 3%, which is typical for such a structure. The point source is positioned at the front of the 
structure. The ray computation is performed using a triangulated surface consisting of 11623 triangles and 
with a 6th order Legendre polynomial basis in direction space. The Nastran grid contains 40 670 elements 
comprising a mix of piecewise linear triangles and quadrilaterals. As would be expected one sees more 
oscillation in the full wave model. The DFM prediction of the overall energy flow in regions of both high and 
low curvature matches well with the Nastran results. The flow of energy along the side walls of the central 
mount, as well as features such as shadow regions due to holes in the structure and channelling effects due 
to variations in the shell thickness can be observed as common in both plots. Such geometrically dependent 
features would be entirely absent from SEA-type models and represent a major advance for high frequency 
vibro-acoustic simulation methods. 

Figure [6] gives the kinetic energy averaged over four different receiver regions, which are displayed as the 
blue regions on the inset structure plots. The response is now shown for a range of damping levels between 
1% and 6%, and for receiver regions with differing levels of proximity to the source point. The correspondence 
between the mean of the 41 Nastran solutions and DFM is remarkably good in Fig[6](a) to (c). The regions 
are spread across the whole structure and provide strong verification of our approach. In Fig.[6](d), the DFM 
results deviate slightly for high damping. Note that this region is close to the region in Fig. [6] (c) and the 
deviations thus reflect local variations in the FEM solution. The results clearly demonstrate that curvature 
related wave-barrier effects are correctly predicted by DFM using local reflection/transmission matrices. 
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The computations in this section were performed in parallel on four cores of a desktop PC in less than 
6 minutes. The code was written in C++ with opcnMP. The large sparse non-symmetric linear system 
of dimension 732 249 which arises has been solved efficiently and accurately using a stabilized bi-conjugate 
gradient iterative solver. 

5 CONCLUSIONS 

We present discrete flow mapping as a highly efficient method for solving stationary phase space flow equations 
in complex domains and apply the method to two illustrative examples. In particular we highlight the 
versatility of the method, its capability to handle large scale models efficiently and emphasise the fact that 
both geometric details and flow directivity are fully included in the calculation at a moderate computational 
cost. In the special case of ray focusing, the method clearly captures the foci albeit not reproducing the 
actual singularity. In the vibro-acoustic example, the full complexity of the model is taken into account via 
the mesh functionality of DFM. Deviations from a pure geodesic flow due to regions of high curvature have 
been reproduced without compromising the efficiency of the method. DFM can thus serve as a practical 
tool for solving phase space transport problems with applications in engineering ranging from acoustics to 
structural mechanics. Extensions of the method to electrical engineering, as well as to fluid mechanics and 
other flow problems are evident. 
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